c
c     ... General Numerical Tools for HOD Calc ...
c

ccc   General tool
   
      real function pois (xlambda,nn) !Poisson Prob Distr.
ccc   Max 17.0
      rnfac=1.0
      if (nn.gt.1) then
         do i = 1,nn
            rnfac=rnfac*i
         end do
      endif
      pois = exp(-xlambda)*xlambda**nn/rnfac
      return
      end

ccc   numerical recipes stuff. 
      SUBROUTINE qromb(func,a,b,ss)
      INTEGER JMAX,JMAXP,K,KM
      REAL*8 a,b,func,ss,EPS
      EXTERNAL func
      PARAMETER (EPS=1e-6, JMAX=26, JMAXP=JMAX+1, K=5, KM=K-1)
!      PARAMETER (EPS=5e-7, JMAX=30, JMAXP=JMAX+1, K=5, KM=K-1)
CU    USES polint,trapzd
      INTEGER j
      REAL*8 dss,h(JMAXP),s(JMAXP)
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd(func,a,b,s(j),j)
        if (j.ge.K) then
           call polint(h(j-KM),s(j-KM),K,0.d0,ss,dss)
           if (abs(dss).le.EPS*abs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
 11   continue
      write (0,'(A,1PE9.2)') 'too many steps in qromb, delta='
     $     ,abs(dss)/abs(ss)
      ss = -1.0
      END
      
      SUBROUTINE trapzd(func,a,b,s,n)
      INTEGER n
      REAL*8 a,b,s,func
      EXTERNAL func
      INTEGER it,j
      REAL*8 del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END

      SUBROUTINE polint(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL*8 dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL*8 den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)
     $         write (0,'(A)') 'failure in polint'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END

      FUNCTION bessj0(x)
      REAL bessj0,x
      REAL ax,xx,z
      DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
     *s1,s2,s3,s4,s5,s6,y
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6
      DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
     *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
     *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
     *651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
     *s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
     *59272.64853d0,267.8532712d0,1.d0/
      if (abs(x).lt.8.)then
        y=x**2
        bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
     *(s4+y*(s5+y*s6)))))
      else
        ax=abs(x)
        z=8./ax
        y=z**2
        xx=ax-.785398164
        bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      return
      END

      SUBROUTINE qtrap(func,a,b,s)
      INTEGER JMAX
      REAL*8 a,b,func,s,EPS
      EXTERNAL func
      PARAMETER (EPS=1.e-2, JMAX=20)
CU    USES trapzd
      INTEGER j
      REAL olds
      olds=-1.e30
      do 11 j=1,JMAX
        call trapzd(func,a,b,s,j)
        if (abs(s-olds).lt.EPS*abs(olds)) return
        olds=s
11    continue
      write (0,'(A,1PE9.2)') 'too many steps in qtrap, eps=',eps
      END
c

c
c
c
      subroutine integ_per (f,xmin,xmax,r,eps0,ans)
      implicit real*8 (a-h,o-z)
      ! Parameters:
      ! f (double precision integrand)
      ! xmax,xmin (min and max bound)
      ! r (2pi/period)
      ! sum_all  (output integral)
      ! The integrand is oscilatory with ~sin(xr) cos(xr) etc...
      ! Integrations are made perid by period, from xmin to xmax
      ! but the end point may not be exactly xmax
      ! but the largest xmin+n*2*pi/r not exceeding b 
      external f
      nst=10
      imax = 12
      eps = 1.0d0
      ans = 9d99
      ans0 = -9d99
      i=0
      do while (eps.gt.eps0.and.i.le.imax)
         ans0 = ans
         call integ_per_nst (f,xmin,xmax,r,nst,ans,0)
         eps = abs(ans0-ans)/abs(ans0) 
         nst = nst*2
         i = i+1
      end do
      if (i.gt.imax.and.eps.gt.eps0)
     $     write (0,'(A,1PE9.2)') 'Warning! eps=',eps
      return
      end
cc
cc   Integration routine specifically designed for sin(xr) etc...
cc
      subroutine integ_per_nst (f,xmin,xmax,r,nst_pi,sum_all,ishow)
      implicit real*8 (a-h,o-z)
      ! Parameters:
      ! f (double precision integrand)
      ! xmax,xmin (min and max bound)
      ! r (2pi/period)
      ! nst_pi (number of steps per pi)
      ! sum_all  (output integral)
      ! ishow (=1, display integration by step to STDERR) 
      ! The integrand is oscilatory with ~sin(xr) cos(xr) etc...
      ! Integrations are made perid by period, from xmin to xmax
      ! but the end point may not be exactly xmax
      ! but the largest xmin+n*2*pi/r not exceeding b 
      external f
      pi = 3.14159265358979323846264d0
      n_per = int((xmax-xmin)*r/2.0d0/pi)
      dx= pi/dble(nst_pi)/r
      sum_all = 0.0d0
      sum_p = 0.0d0
      nstep=0
      if (ishow.ge.1) then
         write (0,'(A,3(1PE11.2))') '!r,xmin,xmax=',r,xmin,xmax
         write (0,'(A,I6)') '!Check r,x0,sum_p,sum_all, n_per=',n_per
      end if
      if (n_per.le.2) then ! Few periods/direct integration
         nstep = max(nst_pi*n_per,nst_pi)
         dx = (xmax-xmin)/dble(nstep)
         do i = 1, nstep
            x=xmin+dx*dble(i)
            sum_all = sum_all+scoeff(i,nstep)*f(x)*dx
         end do
      else
         do ip=0,n_per
            x0 = xmin+ip*2.0d0*pi/r
            if (x0*r.gt.10.0.and.abs(sum_p)*(xmax-x0)/dx*nst_pi*2.0.
     $           le.sum_all*1e-5) then
               go to 3
            else
               sum_p = 0.0d0
               do i = 0, nst_pi
                  x1=x0+dx*i
                  x2=x1+pi/r
                  sum_p = sum_p+scoeff(i,nst_pi)*dx*(f(x1)+f(x2))
                  if (ishow.ge.1.and.ip.eq.0) then
                     write (0,'(5(1X,1PE10.3))') x1,x2,f(x1),f(x2),sum_p
                  end if
               end do
               sum_all = sum_all+sum_p
                                ! Check
               if (ishow.ge.1) then
                  write (0,'(4(1X,1PE10.3))') r,x0,sum_p,sum_all
               end if
            end if
         end do
      end if
 3    continue
      if (ishow.ge.1) then
         write (0,'(A,I8,1PE11.2)') '!Check...end at ip,x=',ip,x0
      end if
      return
      end
c     
c
      real*8 function scoeff (i,nmax) !__ coeffs for Simpson int.
c     Give coefficinets for the Simpson integration
c      0 <= i <= nmax   
      ii=min(i,nmax-i)
      if (nmax.ge.8) then
         if (ii.eq.0) then
            scoeff=17./48.
         else if (ii.eq.1) then
            scoeff=59./48.
         else if (ii.eq.2) then
            scoeff=43./48.
         else if (ii.eq.3) then
            scoeff=49./48.
         else
            scoeff = 1.
         end if
      else if (nmax.ge.4) then
         if (ii.eq.0) then
            scoeff=5./12.
         else if (ii.eq.1) then
            scoeff=13./12.
         end if
      end if
      return
      end
c
      subroutine intra_mono (tbl,val,imin,imax,ff)
!     Given a real array tbl of monotonically increasing/decreasing order,
!     return imin,imax s.t. val is between tbl(imin) and tbl(imin+1), 
!     imax=imin+1
!     ff is ((val-tbl(imin))/(tbl(imax)-tbl(imin))
!     Input imin, imax defines the region where the table is defined.
!     1/2 folding method. If val is outside of the range,
!     ff<0 or ff>1, giving extraploation from the last two values.
      implicit none
      real tbl(*)
      real val, ff, pm
      integer i, imin,imax
      pm = sign(1.0,tbl(imax)-tbl(imin)) ! Increasing or Decreasing
      do while (imax-imin.gt.1)
         i = nint((imin+imax)/2.0)
         if (pm*tbl(i).gt.pm*val) then
            imax = i
         else
            imin = i
         end if
      end do
      ff = (val-tbl(imin))/(tbl(imax)-tbl(imin))
      return
      end 
c
      subroutine intra_mono_d (tbl,val,imin,imax,ff)
!     Given a real*8 array tbl of monotonically increasing/decreasing order,
!     return imin,imax s.t.  val (real*8) is between tbl(imin) and tbl(imin+1), 
!     imax=imin+1
!     ff (real) is ((val-tbl(imin))/(tbl(imax)-tbl(imin))
!     Input imin, imax defines the region where the table is defined.
!     1/2 folding method. If val is outside of the range,
!     ff<0 or ff>1, giving extraploation from the last two values.
      implicit none
      real*8 tbl(*), val
      real ff
      integer imin,imax
      integer i
      real pm
      pm = sign(1.0d0,tbl(imax)-tbl(imin)) ! Increasing or Decreasing
      do while (imax-imin.gt.1)
         i = nint((imin+imax)/2.0)
         if (pm*tbl(i).gt.pm*val) then
            imax = i
         else
            imin = i
         end if
      end do
      ff = (val-tbl(imin))/(tbl(imax)-tbl(imin))
      return
      end 
c
      character*(*) function cut_space(rline)
c     Removing front spaces
      integer nmax
      character*(*) rline
      character*1 tab
      tab=char(9)
      nmax=len(rline)
      cut_space=rline
      i=0
      do while ((cut_space(1:1).eq.' '.or.cut_space(1:1).eq.tab)
     $     .and.i.lt.nmax)
         cut_space=cut_space(2:)
         i=i+1
      end do
      return
      end
c
      subroutine spline(x,y,n,yp1,ypn,y2)
      integer n,NMAX
      real*8 yp1,ypn,x(n),y(n),y2(n)
      parameter (NMAX=500)
      integer i,k
      real*8 p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END


      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL*8 x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL*8 a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) then
         write (0,*) ' !! *** Error in splint'
         stop
      end if
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END





